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Abstract. It is shown that the density of the ratio of two random variables with the same 
variance and joint Gaussian density satisfies a non stationary diffusion equation. Implications of 
this result for kernel density estimation of the condensed density of the generalized eigenvalues of a 
random matrix pencil useful for the numerical inversion of the Laplace transform is discussed. 
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Introduction. The density of the ratio of two random variables with joint bivari- 
ate Gaussian density has been derived by several authors and it is important in many 
applications (see e.g. El [HI US]). In the sequel it is proved that, when the two 
variables have the same variance, this density satisfies a parabolic partial differential 
equation whose coefficients depend on both the independent variables. The proof is 
based on standard properties of the confluent hypergeometric functions of the first 
kind. A motivation for deriving such a PDE is provided by the problem of the numer- 
ical inversion of Laplace transform from noisy discrete data [U [3] . This is a classical 
ill-posed problem. Insights for its stable solution can be obtained from knowledge 
of the marginal densities of the damping factors of a multiexponential model which 
represents a discretization of the Laplace transform. This problem can be restated 
in terms of the condensed density of the generalized eigenvalues of a matrix pencil 
built from the observations. In a recent paper [5] an adaptive kernel density estimator 
based on linear diffusion processes has been proposed which has several advantages 
over the existing methods. In the sequel a kernel density estimator in the class con- 
sidered in [5] , based on the proposed diffusion equation, for estimating the condensed 
density mentioned above is proposed. A Montecarlo simulation allows to appreciate 
its merits with respect to a Gaussian kernel estimator and its effectiveness for the 
numerical inversion of the Laplace transform. 

The paper is organized as follows. In the first section the density of the ratio of 
two random variables with joint bivariate Gaussian density is shortly derived in terms 
of confluent hypergeometric functions of the first kind. In the second section the 
PDE is derived. In the third section the kernel density estimator based on the PDE 
is derived and the conditions which need to be met by the function whose Laplace 
transform has to be inverted in order to get good results are specified. In the last 
section the merits of the proposed method are shown by a MonteCarlo simulation. 

1. The density of the ratio of two jointly distributed Gaussian variables. 

Let us assume that the random variables (v, w) have a joint Gaussian density 
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with |S| > and the mean, the covariance matrix and its inverse are given by: 
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Let \Fi[a,j3, z] be the confluent hypergeometric function of the first kind. The fol- 
lowing lemma holds: 

Lemma 1.1. If a £ R + , beM,Vn€N 
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if n is odd 
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But (see e.g. [i 3.462,1]) 

X n+1 g(X)dX = (2a)- 2 >-r(n + 2)e^ J D_ ( „ +2) 
where the parabolic cylinder function Z3_( Jl+2 )(z) is given by 
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We can get the density of the ratio of Gaussian variables as a simple consequence of 
this Lemma (see also [TB]): 

Theorem 1.2. 7/E > ; i/ie density of the ratio x = — is given by: 
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and a(x) > Vx, c > 0. 

Proof. The density of the ratio x = — can be written as: 
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By the change of variables A = u, A* = ^ with Jacobian |A| we get 
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Moreover a(x) > Vx because the quadratic equation in x 

a 2 w - 2 7 x + a 2 v x 2 = 
has no real roots as E > 0, hence, by Lemma H. 11 with n = 
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Finally we notice that c > as cr 2 f 2 — 2^v w v v + <r 2 . l 'w > [p w v v — (t v u w ) 2 > 0, because 
E > 0. □ 

Corollary 1.3. If v v = v w = ant! E = cr 2 / i/ien /z(a;) = ^ct^ . 
Proof. In the considered case we have 
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, b(x) = 0, c = 0, |E| = er 4 , 
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2. The diffusion equation for the density of the ratio of two jointly 
distributed Gaussian variables. Let us assume that 
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and define t = a 2 . By making explicit the dependence on t in a(x),b(x),c, \H\,h(x) 
we get 
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Therefore if v v ^ and a = — we have 
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We have 

Theorem 2.1. 
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in the weak sense. 

Proof. The density h(x,t) can be rewritten as 
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Taking the limit t — > oo in this expression we get the first equality in the first part of 
the thesis. The second equality is obtained by substituting v v = v w = in equation 
(|2.2p . To prove the second part we notice that 



e 2t ( l2 - 2l " +1 ) v v (l-px) + v w (x-p) 
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holds for all continuous compactly supported functions -F, and so t) converges 

weakly to 5 (x — in the sense of measures ( [TSJ Theorem 1.18]). □ 

The properties of h(x, t) stated above suggest, when v v ^ 0, the existence of a diffusion 
equation ruling the behavior of h(x, t) for varying t (when v v = v w = 0, h(x, t) does 
not depend on t). To prove that this is indeed the case we need the following Lemmas: 
Lemma 2.2. If v v 7^ and \p\ < 1 then 
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By the change of variables A = u, /i = ^ with Jacobian |A| we get 
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where /(A) = e -a(=M)A 2 +26(x,t)A-c(t) _ e -c(t) g ^ and a (a;,t),6(a;,t),c(i) are given in 
equations (|2.1j) . and 
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In the same way we get 

h x (x,t) = — \A^(x,t) j \X\X 2 f(X)dX + B^(x,t) J \X\Xf(X)dX 
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By using the same notations of Lemma ITTTI we get the thesis. □ 
Lemma 2.3. If a e M+, b e 2R, 

L 3 (M) = Wi£i(a:,t) + W 2 L 2 (x,t) 
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and 



where 
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Proof. By Lemma 11.11 we have 
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By [TJ 13.4.3] we have: 
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We can now prove the main theorem: 

Theorem 2.4. If v v 7^ and \p\ < 1, the density h(x,t) solves the partial 
differential equation 



(2.4) 
(2.5) 



h t (x, t) = T>h(x, t) 
V = 



e« t * t )i +S | t ), 



dx 



where the diffusion coefficient is 
(2.6) D(x,t) 
the source coefficient is 



and the convection coefficient is 
P 1 (x)+tP 2 (x) 



C(x,t) 
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P^x) P 3 (x)(Q[(x) + Q' 2 (x)) 



t(Qi(x) +tQ 2 (x)) Q!{x)+tQ 2 (x) (Qi(x) +tQ 2 (x)) 2 



where 



Pi(x) = 2(u w ~ v v xf\v v + v w x - (y w + v v x)p](p 2 - 1) 
P 2 (x) = (1 + x 2 - 2xp)[v w {p- x)(3x 2 -6xp+ Up 2 - 8) + 

v v (2 - 9x 2 + lOxp + 3x 3 p - hp 2 - xp 3 )] 
P 3 (x) = (1 + x 2 - 2xp) 2 {v w {\ -x 2 + 2xp - 2p 2 ) + u v [p + x(-2 + xp)}} 
Qx{x) = 2(1 - p 2 )(v w - v v x) 2 (v w - v v p) 

Q 2 (x) = 2(p 2 - l){v w (l + Ax 2 - 8xp + 3p 2 ) - u v [p + x(3x 2 - 5xp + p 2 )}}. 

Qi(x) +tQ 2 (x) is a cubic polynomial with one, two or three real zeros depending on 
the values of t, v v ,v w ,p. 

Proof. Dropping the dependencies on (x, i), by Lemma 12.21 we have 
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By Lemma \2. 2 1 we have 
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where 



(j(xx) _ A (xx) + p(xx) w ^ + ^Os)^ £)(:ee) _ F^Wl + E^Wz 



We get 



L 2 = e c(t) 27r 
Li = r''"'2;r 



A(x) D (xx) _ B (x)(j(xx) 



Substituting these expression in 
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and remembering that 



we get 
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where 
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Substituting the expressions for A'*' , A' 1 ' , , B W , C<"), D' 11 ' , given 
in Lemma 12.21 and noticing that 

dG 

C(x,t) = G x and D(x,t) = G xx 

ox 

we get the expressions reported above. Moreover Qi(x) + tQ 2 {x) — is a cubic 
polynomial equation whose discriminant can be positive, negative or zero depending 
on the values of t, v v ,v w ,p. □ 

3. A density estimation problem. Let 

Lf(s) = / f(t)e- st dt 
Jo 

be the Laplace transform of a function f(t) € Li(lR + ). Let us denote random quan- 
tities by bold characters. Let be 

d k = L f (kA s ) + e k , A s >0, k = l,...,n 

where e k are i.i.d. Gaussian zero mean random variables with variance er 2 and let us 
consider the problem of making inference on f(t) from R independent realizations of 
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d = [di, . . . , d„]. The problem can be severely ill-posed. An approach to its solution 
consists in approximating the Laplace transform by a finite sum, assuming n even 



L f (kA s ) » fj e ~ aj{k ~ 1} > a 3 > °> P 



3=1 

and in solving for the unknowns {fj, otj}, j = 1, . . . ,p in the multiexponcntial model 
(for simplicity the same symbols are used): 

<U = £ fie- **- 1 ) +e k = J2 tet 1 + "■■ 

3=1 3=1 

In the noiseless case the problem consists in interpolating the data 

v 

(3.1) Sk = Y J h^ a ^ 1 \ k = l,... >n 

3 = 1 

by means of a linear combination of real exponential functions Q(t) = e~ ajt , j = 
1, . . . ,p. To this aim let us consider the Hankel matrices 

U (s) = U(s , . . . , s n _ 2 ), Ui(s) = U(si,...,s n -i) 

where 



U(xi,. . . ,x„_i) = 



X2 X 3 
Xp Xp-\-\ 



■'■ p 

Xp+1 
Xn—1 



It is well known (e.g. [7]) that, provided that det([/o) ^ 0,det([/i) ^ 0, a unique 
solution exists. If £ and W denote the generalized eigenvalues and eigenvectors of the 
matrices (Ui, Uq) then the solution is given by 

c=i, i=w t s = v{q- i s 

where V(£) is the square Vandermonde matrix based on £ and T denotes transposi- 
tion. Hence the critical quantities which the solution depend on are the generalized 
eigenvalues £. They can be computed by the generalized Schur decomposition of the 
matrices {U\,Uq) [5]: 

Ui = QSZ T , U = QTZ T 

where Q and Z are orthogonal matrices, and S and T are upper triangular matrices 
such that £j = 5"-. In the noisy case the matrices Uo, Ui are random and the gener- 
alized eigenvalues £j, j = 1, . . . ,p are random variables. Their marginal densities are 
all equal to the their condensed density (see e.g. [2j Lemma 2.4]) which is defined as 



(3.2) 



H{x) = E 



-E^-^3-) 

P ~^ 



3 = 1 



Knowledge of the condensed density is therefore of main importance for making in- 
ference on the generalized eigenvalues £. 
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In a more general context this problem was studied in 3 a where a stochastic 
perturbation method for estimating the condensed density (13.21) based on a single 
realization of d was proposed. Here we assume to have R independent realizations 
S r \ r = 1, . . . , R of d and we are seeking a kernel estimator of the marginal densi- 
ties. In a recent paper [5] it has been shown that kernel estimators based on parabolic 
partial differential equations can be considered and the underlying PDE can be used 
to estimate the optimal bandwidth and to take into account some kinds of prior infor- 
mation through suitable boundary conditions. Gaussian kernels belong to this class 
as they satisfy the heat equation. In the specific case considered here the Gaussian 
kernel estimator of (|3.2I) takes the form 

(3-3) A G (M) = ^£-]-i>(z4 r) ^) 



where 



r=l ^ r k=l 



where Q , k = 1, . . . ,p r are the real generalized eigenvalues of {u[ r \ Uq ) built from 
(discarding the complex conjugate pairs). It turns out that Hq{x, t) is the unique 
solution of the diffusion equation 



with initial condition Hq(x, 0) = H e (x) where 

lx i Pr , i 

r=l k=l 

is the empirical condensed density of the generalized eigenvalues. 

We now notice that if p = 1 the only generalized eigenvalue £ = d 2 /d! is the 
ratio of two uncorrelated Gaussian random variables with the same variance t = a 2 
and mean fxC^ and /i respectively and its density was derived in Section 1. Moreover 
in Section 2 a diffusion equation was derived which is satisfied by this density. The 
idea is then to replace the standard diffusion operator which gives rise to a Gaussian 
kernel density estimation with a more specific diffusion operator related to the one 
defined in Theorem [231 However we can not use straightforwardly the operator (|2.5|) 
because the theory developed in [5] holds for diffusion operators with coefficients 
independent of t and positive diffusion coefficient. On the other hand when p > 1 the 
generalized eigenvalues are the ratio of variables which are not Gaussian. Therefore 
in any case, when proposing a modified operator based on (|2.5p . we are looking for 
a suboptimal solution to the kernel selection problem. However it turns out that the 
generalized eigenvalues can be approximated by the ratio of Gaussian variables and 
the approximation errors of the numerator and denominator are random variables 
whose expectation and standard deviation are proportional to 
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This will be proved in Theorem 13.31 Hence the approximation can be very good if 

rr im 

the signal-to-noise ratio, measured by ~i , is large enough with respect to the 

relative distance of the numbers £i, i — 1, . . . ,p, measured by rii<j(Ci — (j) 6 - 

A modified operator can be built as follows. We first notice that the difficulty 
of the Laplace inversion problem strongly depends on the relative position of the 
Q, j = 1, ... ,p which the interpolation of the noiseless data is based on. Simplistically 
the closer they are the worse the conditioning of the problem is. We then prove that 
in these difficult cases the diffusion coefficient of the operator (|2.5I) is positive in a 
neighbor of the interesting region of the density for a small enough. This is proved 
in Theorem 13.21 We first need the following 

Lemma 3.1. The generalized eigenvalues of the random pencil (Ui,Uo) built from 
the data d— [di, . . . , d„] are given by Do = diag(t_), T>± = diag(t_) ■ diag(C) where 

£j,fj,j = l,...,p are random variables such that dk(uj) — Xw=i /i( w )Cj' (^0) & = 
1, ... ,n, Vw € Q, where O is the space of events. If the generalized Schur decomposition 
of (Ui,U ) is given by 

Ui = QSZ T , U = QTZ T 

then 

diag(S) = Di, diag(T) = D . 



Proof. Let fj,Cji j = 1, ... ,p be the solution of the exponential interpolation 
problem which exists and it is unique a.s. because dct(U ) ^ a.s and dct(Ui) ^ 
a.s. [IT]. If V is the Vandermonde matrix Vjj = C,\ then (see e.g. [3]) 

U = VD V T , Ui = VD!V T . 

But then 

UiV" T D = U V" T Di. 

Therefore the pairs ({jCj,fj),j = l,.,.,p are representatives of the projective form 
[T5] of the generalized eigenvalues of (Ui, Uo) and the thesis follows. □ 
Theorem 3.2. If Pj — corr(Sjj,Tjj) then for h^k 

lim lim oa = 1, Vj 

ICh-Cfc|->Oc7^o rj 

and it exists an open interval I C IR + such that ^ £ I and D(x,t) > 0, x G /, Vt. 
Proof. By Lemma 13.11 we have 



Pi 



(3.4) 



E 





E[S n ] 


T - 


E[T n ] 








-E[Tn}) 2 } 






f,- 


m 


Jmc 3 


-mc 3 }) 2 } 




-E[fj]) 2 ] 



For each realization, and fj- are analytic functions of d^ in a small neighbor of s 
( [3j Lemma 2] ) , therefore they admit Taylor series expansions around s 



n ^ n 

= + X] 9i€i + o E C ih e,e h + ... 



i=l 



i,h=l 
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n 1 n 

f j = u + X c%ei + o X G J i/ie 



2 / ^ " jcif*~c^h 
i=l i,/i=l 

Truncating after the first order terms and substituting these expressions in (|3.4j) . after 
some long but simple calculations, we get 

15 _ „ ELi c ? + /j c 'g' 

° ? w /v * — ^ W. O n -i — ^ 7). \ A « O / i — ^ T7. 1 i — ^ T7. O / 1 — ^ TJ. \ Z \ ' 



{(0 ELi c ? + /j ELi c i9i? + ff (e?=i c ? E?=i - (ELi c i5i) 

But if for some h ^ k, \Ch ~ Ck\ Q then |c,| — >• oo Vi. In fact 



9di | d=£ ddi [d=2 

and 

0(e?Y(C)- r s) a(ELi^(0« fc ) ^ ^(0 ^ 

Ci = — ^ — = a^— -g^r sfe+ ^) 

where Ujfc = eJV~ T e fc . But (see e.g. [13]) 
and therefore 

Ci = - X) ffc^J r ) - fe + ^ ( -' 

The first part of the thesis then follows by noticing that 

detF= Xl(a-Cfe) 



and 



jj Ei=l C i + fj Ei=l C i9i _ , 



(|ci|,...,|c„|)->oo 



{(0 Er =1 4 + h EF =1 ^) 2 + /| (Eti EILi 3? - (EIU C ^) 2 ) } 

because £j > 0. To prove the second part, let us consider the Taylor first order 
approximation of the diffusion coefficient around p = 1 and x = i ^ L : 

D(x, t) = A(t) + ( x -^-) (B(t) + 0(1- P f) +0(l-p) 2 + o(x-'^) 2 
where 

(y v - v w f (u v - v w f {yl + 6v v v w + v^) (1 - p)[y v + v w f 



A® 



A(l-p)(ti4) Mv$ 16 (t«*) 
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and 



B(t) 



7(v v - Vw) 3 (y v - v w ) {yl + 28VwUy + Ivl) 



1) 



8(^3) 



l&tiA - Votvli 



But A(t) > as v v , v w have the same sign because Q > 0. Therefore we get the thesis 
by the permanence of sign theorem. □ 

By using Theorem 13.21 we can define the modified operator as the operator (|2.5[) 
where the coefficients are evaluated at a fixed suitable value to. When p = 1 the 
variable t represents the common variance of the numerator and denominator of the 
generalized eigenvalue. In order to choose to we can then look for the element in the 
set of densities (|2.3[) which best fits the empirical condensed density H e (x), i.e. 

(*oA) = a,rgmm te \\h(x,t;9) - H e {x)\\l 

where 9 — {^Sp}- Let us denote by 



■Do = 



1 1 0)9 * l - + C (x,t )— + 5 to • 
ox 



dx 



this modified operator and define the kernel estimator 
(3.5) 



H P {x,t*) = ^Y,±-Y,h^{x,t*) 



R 1 — ' p r 

r=l l r k=l 



where 



• h( r ' k ^ [x, t*) is obtained by equation (|2.3j) by replacing ^ by obtained 

V ^kk 

by computing the generalized eigenvalues by the Schur decomposition of the 
matrices (ll[ r \ Uq^) built from d^, taking the p r real ones (discarding the 
complex conjugate pairs), and by replacing p with the sample correlation 



(r) (r) 

coefficient p of the pooled real S kk ' , T kk , r = 
the optimal bandwidth is given by 5, eq.23] 



1 />': fc = l,-. 



• Pr 



E 



(yD(x,toj) 



\ 



2/5 



where E 



2R^\\V h(x,t )\\l 
is estimated by 



R * — ' p r 

»=1 ^ r k=l 



\T>oh(x,to)\\2 is estimated by 



R £ — ' Pr 

i=l rr k=l 



and \\h\ r ' (x, ^o) || 2 is computed by numerical quadrature; 
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• Z)( r ' fe ) [x, to) denotes the diffusion coefficient computed by replacing in formula 
(|2.6[) v w ,is v ,p by the the same values used for hS r ' k ^ (x, t*). With the same 
substitutions we obtain h[ r ' k ^ (x, to) by formula (|2.7[1 : 

By the second part of Theorem 12.11 Hp(x,t) is the unique solution of the diffusion 

equation 

(3.6) ^- t Hp(x,t)=V H P (x,t) 

with initial condition Hp(x,0) = H e (x). 

In the next Theorem conditions under which the distribution of the generalized 
eigenvalues is well approximated by the distribution of the ratio of Gaussian variables 
are specified. 

Theorem 3.3. The generalized eigenvalues (ijCj,ij),j = l,...,p o/(Ux,Uo) 

are given by 

n 
i=l 



f^fj+^cj^ + vP (x) 



where hji and Cji do not depend on f,xisa point of M n lying in the interior of the 
line segment joining d and s, and 

^(s)} < , NB , h=l,2 



2 IlUlfrUr^r-Q 6 



var[ ^ ( - )] - Tni./m^fc-o 12 ' 

where Fh(-), h = 1, 2 are polynomials in . 

Proof. Let $ : J?" — » 2R" be the map that associates to each n— vector the n/2 
pairs corresponding to the projective form of the generalized eigenvalues of the pencil 
(Ui, Uq) built from the n— vector. It was proved in [21 Lemma 2] that <!> is analytic. 
We can then consider the first order Taylor series expansions with remainder of Cj 
and fj, as functions of d, around s ([17j Th. B]): 



n 1 71 



2 <H 

i — 1 — 1 



n 1 n 



2 

i,h=l 



where 
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ddi | d=s ' 11,h ddidd h i d=a ." 



We notice that g 3i = ^ and analogously for c,-j , , G^ . Let us denote by / 
and #> = ^.Lt be 



(h) 

ji 





' 1 


1 


.. 1 




Ci 




■• C P 


V = 


C? 


CI • 


.. c 2 

Sp 




1 

. SI 


/-n— 1 
S2 


/-n—1 
■ • Sp 



and 



V, 



(i) 



At) 

Sit 

2 CiCh 



/•(I) 
S2i 

2 C2C 2l 



(»-i)cr- a ci? (^-i)c 2 "- 2 c^ 



2GC 



(i) 

■pi 



By derivating both members of equation ()3.1|) with respect to Sj we have 



(3.7) 



a§. 

8S: 



= e„- 



ds, 



vPf + VfW^VDgf + VfW 



where D^) is the diagonal matrix built from Cu > 



V 




1 

2Ci 




1 

2C2 



(n-l)Cr 2 (n-l)^ 



■•-Cff and 




1 

2C P 

(»-l)CS 



But then if 



o = [f,C T } T 

and D/ is the diagonal matrix built from fi, ■ ■ • , f p , equation (|3.7[) becomes 



(3.8) 
and therefore 

We then have 



- WOP = e. 



Cji = ej [I \0}W~ 1 e l , 9ji = ej [0 : I] W~ 1 e t . 



dsi 



eJ(C ] [l':0} + f 3 [0\l])W- 1 e i . 
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But 



hence 



W=[V\V] 



Cji =ej[l\0] 



' I 







' / 








D f . 


, w- 1 = 








/ 
D 



f 



[V-.Vr'e^ejil-.OW-.V]- 1 ^ 



is a function of £ only (it does not depend on /). As 

h j i = ([Qej\0} + [0\ej])[V\V}- 1 e i 

does not depend on / we get the first part of the thesis. 
Let be 



Gjih — 



d% 



r< _ u h 

l-'jih — 



ddidd h | d=£ 

where, for simplicity, the same symbols as before were used, and let be 



1 ™ 



,/i=l 



where 



because 



But then 



(2) 



1 ™ 



i,h=l 



Mid 



dd 



i |d=s 



S^^U)] = y tr(^), £fof |s] = y tr(Q) 



and, by Isserlis's theorem, 



4 



2 
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4 ? ' 



■2 

jih 



To conclude the proof we need an expression for Cjih and Gjih ■ If 

d 2 e 



ih 



dsidsh 



i/iJ 



by e.g. [13l Ch.5] we have 

Can = ej[I : 0}T ih = ej[I : Q] A^-ig.] = _gT[/ | o]^- 1 ^^- 1 ^ 



and 



where 



and 



G 



w v = d -K = 

ds h 



ds h 



VD$\VD$ +VD S D$ 



V = 



(n-l)(n-2)Cr 3 (n-l)(»-2X5 





2 

(n-l)(n-2)C 



n— 3 
P 



and and £>i^ are the diagonal matrices built respectively from Cu ■> ■ ■ ■ i Cpi^ an d 



/ii > • • • j /pi J - We now notice that the elements of [V : V] 1 are rational functions of 
Ci, • . • , Cp- More specifically by [TO] 



f (i) 



[V : V]~ l = Dw ■ X, D w = 



D 3 
D 2 



where 



D = diag 



life - CO, 11^ -C P ) 

i^l i^p 



and the elements of X are polynomials in £i> . . . , Cp. But 



diag 



jC 1 ) 



= diag [W^eJ = 
I 







(i) 



D 



/ 



'C* 

diag [XeJ 



Z)- 3 
D- 2 

D- 3 D xft 

D- X D- 2 D XQ% 
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and therefore 



VD^D- 2 D xii \VD- 3 D xfi + VD~ 2 D XCi 



VD x(i \VD xfi + VDD x(i 



D- X D- 2 

D~ 3 







Let us consider the matrix equation in the unknown B 

WB = W {1) . 

As the right block of W and the left block of are both equal to V times a 

diagonal matrix, B must have the form 



and 



B 



B\2 

B22 



B 12 
Db B22 



D B = Dj 2 D- 2 D XCt 



W 



1 (vD- 3 D xft + VD- 2 D x( ^ 



I 
Df 



\y':V]~ 1 VD~ z D X f i 



I 





D 



f 



D- 3 D Xfl + 



D- 3 

D } lD ~ 2 



[V:V]~ 1 VD~ 2 D x ^ i 
XD- l D XCl 



and XV = XD 



because it turns out that 

[v'-.v^v = 

and the elements of X are polynomials in Ci, . . . , C P ; therefore 



B 12 = D~ 3 X\D~ 1 D X (-i 



and 



But then 



B22 = D^D- 3 D xfl + DJ 1 D- 2 X 2 D- 1 D XCi 



W-^W- 1 = BW- 1 = 
os h 



B 12 
Db B 2 2 



D- 3 

D j lD ~ 2 



Xn X12 
X21 X 2 2 



B^D^D- 2 
D B D- 3 B 22 D- 1 D- 2 





' X n 


X\2 




X21 


X22 



21 



B 12 DJ L D~ 2 X 21 B 12 Dj L D- X 22 

D B D- 3 X n + B 22 Dj 1 D- 2 X 21 D B D- 3 X 12 + B 22 D^D- 2 X 22 



U n U 12 

U 2 \ u 22 

where 

U n = D- s X 1 D- l D XCi DfD- 2 X 21 
U 12 = D^X.D-'Dx^Dj'D^X^ 

U21 = D B D- 3 X n + (Dj 1 D- 3 D X f i + Dj 1 D- 2 X 2 D- 1 Dx C i) D^D~ 2 X 2l 

C/22 = D B D- 3 X 12 + (p^D- 3 D Xfl + Dj 1 D- 2 X 2 D- 1 Dxci) BfD~ 2 X 22 

Remembering that Ci> ■ • ■ > Cp e (0, 1), it follows that Cj%h an d Hjih are rational func- 
tions such that the numerators are polynomials in Ci > • ■ • > Cp an d a lower bound for 
the denominators is flr=i f r IIr^s(Cr ~ Cs) 6 because some further simplification of 
common factors such as (£ r — £ s ) in the numerator and denominator can occur. This 
fact follows easily for Cjih while for Hjih we remember that 

Hjih = ^jihCj ^Cji9ji ~t~ fj^jih 

and we notice that U 2 \ and U 22 are left multiplied by DJ 1 , therefore fjGjih is a 
rational function such that the numerator is a polynomial in Ci > ■ • ■ 1 Cp an d a lower 
bound for the denominator is J\r=i f r Wr^s^ r ~ Cs) 5 - Moreover as Cji,gji do not 
depend on / the claim follows as well as the thesis. □ 

As a final remark we notice that when the densities of f,-£ - and fj are approx- 
imately Gaussian, also their joint density is approximately Gaussian, because the 
density of £ -|f,- is approximately Gaussian too. 

4. Simulation results. In order to illustrate the possible advantages of the 
proposed kernel estimator, the following MonteCarlo simulation was performed. N = 
10 5 independent realizations of a noisy multiexponential signal of length n — 126 with 
three components 

3 

d fc r) =XX- -1 +eL r) , C= [0.8,0.9,0.95], cr = 1.5-10- 3 , k = 1, . . . , n, r = 1, . . . , N 

3=1 

were considered. For r = 1, . . . , N the p — n/2 generalized eigenvalues were com- 
puted as well as their empirical condensed density that was taken as the refer- 
ence distribution that we want to estimate starting from the first R — 250 samples 

(r) 

d k , r = 1, . . . , R. The noise standard deviation a was chosen large enough to make 
at least one of the three modes hardly detectable by visual inspection in the em- 
pirical condensed density based on R = 250 samples and small enough to make the 
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three modes visually detectable in the reference condensed density. The number of 
observations was chosen as a function of a by the rule 

n — argmin fe {fc| \d^\ < <r} 

as a compromise between the opposite requirements of a large sample size and a small 
total noise. 

In the top part of Fig.l the reference distribution evaluated in 256 bins of equal 
size in the interval (0.75, 1) was plotted (right) as well as the empirical condensed 
density based on the first R — 250 samples (left). The kernel estimator Hg(x, t + ) was 
evaluated in 256 equispaced points in the interval (0.75, 1) where t + is the estimated 
optimal bandwidth; the software downloadable by [3D] was used and the result is 
plotted in Fig.l (bottom left). The kernel estimator p.5|) was evaluated in the same 
points and plotted in Fig.l (bottom right). The estimated bandwidths were to = 
1.1 ■ 10 _1 , t + = 1.12 • 10~ 2 . We stress that in this problem what matters are the 
modes of the density because they are estimates of the generalized eigenvalues. A 
smooth estimate with the correct number of modes even if slightly displaced w.r. to 
the true values is much better than an estimate with many modes not related to the 
true ones. Therefore we can conclude that the proposed estimate is much closer in 
a suitable Sobolev norm to the reference distribution than that based on standard 
diffusion. Moreover if we compute the relative maxima of the proposed estimate 
above e.g. a threshold r = 2 we get the modes [0.82, 0.88, 0.95] which are reasonable 
estimates of the true values £ = [0.8,0.9,0.95]. 

To stress the proposed method, a second example was considered where the signal 
has more and closer components. Moreover a was chosen large enough to make one of 
the modes visually undetectable even in the reference density. The multiexponential 
signal of length n — 324 with five components was considered: 

C = [0.88, 0.9, 0.91, 0.92, 0.94], / = [1, 10, 10, 10, 1] 

(7 = 2- 10~ 9 , k = 1, . . . , n, r = 1, . . . , N 

As before, R — 250 samples were used. All the distributions were now evaluated in 
2 13 points in the interval (0.85, 0.96) and plotted in Fig. 2. The estimated bandwidths 
were to — 3.8 ■ 10~ 3 , t + = 1.75 • 10~ 4 . In the reference distribution one mode is lost, 
while the relative maxima above the threshold r = 2 of the proposed estimate are 
[0.880, 0.902, 0.907, 0.919, 0.940]. 

5. Conclusions. The mathematical structure of the density of the ratio of Gaus- 
sian variables given by a partial differential equation has been revealed and exploited 
to solve a classical ill posed problem. The quality of the solution is definitely better 
than the one provided by classical methods. Moreover it turns out that, given a sam- 
ple of observations of moderate size, the quality of the solution can be better than the 
one obtained by a very large sample. The results are apparently robust with respect 
to the Normality hypothesis. It is reasonable to expect that similar benefits can be 
obtained by exploiting the mathematical structure for solving other problems where 
the ratio of random variables plays an important role. 
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Fig. 5.1. Model 1. Top left: empirical distribution of the generalized eigenvalues based on 250 
data samples, evaluated on 256 bins of equal size; top right: empirical distribution of the generalized 
eigenvalues based on 10 5 data samples evaluated on the same bins; bottom left: Gaussian kernel 
density estimation based on 250 samples, evaluated in 256 equispaced points; bottom right: proposed 
kernel density estimation based on 250 samples, evaluated in the same points. 



25 





0.86 0. 



0.9 0.92 0.94 0.96 




0.88 0.9 0.92 0.94 0.96 



Fig. 5.2. Model 2. Top left: empirical distribution of the generalized eigenvalues based on 250 
data samples, evaluated on 2 13 bins of equal size; top right: empirical distribution of the generalized 
eigenvalues based on 10 5 data samples evaluated on the same bins; bottom left: Gaussian kernel 
density estimation based on 250 samples, evaluated in 2 13 equispaced points; bottom right: proposed 
kernel density estimation based on 250 samples, evaluated in the same points. 



